Etudiants : Doré Maxime, Pillard Cécile, Bennis Zineb, Debbagh Boutarbouch Yasmine, Jacquiet Morgane
defaultW <- getOption("warn")
options(warn = -1)
library(lhs)
library(sensitivity)
library(plotly)
library(DiceKriging)
library(randtoolbox)
library(DiceDesign)
options(warn = defaultW)
#Cette fonction prend en argument une matrice X qui contient 3 colonnes de longueur n
# des valeurs réparties uniformément entre -pi et pi
#Les colonnes correspondent respectivement à x1,x2 et x3
ishigami <- function(X, coeff){
return(sin(X[,1]) + coeff[1]*sin(X[,2])^2 + coeff[2]*X[,3]^4*sin(X[,1]))
}
#Cette fonction calcule les indices de Sobol Analytiques
# du premier et du second ordre pour la fonction Ishigami
ishigami.sobol.indices <- function(coeff){
V = 1/2 + coeff[1]^2/8 + ((coeff[2]^2)*pi^8)/18 + (coeff[2]*pi^4)/5
V1 = (1/2)*(1+coeff[2]*((pi^4)/5))^2
V2 = (coeff[1]^2)/8
V3 = 0
V12 = 0
V13 = coeff[2]^2*(pi^8)*(8/225)
V23 = 0
# Initialisation des coeffs
S1 <- St <- array(0,3)
# Coeffs premier ordre
S1[1] = V1/V
S1[2] = V2/V
S1[3] = V3/V #=0
# Coeffs second ordre
S12 = V12/V #=0
S13 = V13/V
S23 = V23/V #=0
S2 <- c(S12, S23, S13)
# Coeffs totaux
St[1] = (V1+V13)/V
St[2] = V2/V
St[3] = V13/V
list(Coeff_Ordre1=S1, Coeff_Ordre2=S2, Coeff_Ordre_Tot=St)
}
n = 4000
X = matrix(NA, n, 3) #Creation d'une matrice vide
x1 = runif(n,-pi,pi)
x2 = runif(n,-pi,pi)
x3 = runif(n,-pi,pi)
X[,1] = x1
X[,2] = x2
X[,3] = x3
#initialisation des coefficients A et B de la fonciton d'Ishigami:
coeff = c(7, 0.1)
cat("Indices de Sobol pour la fonction Ishigami : \n")
(ishigami.sobol.indices(coeff))
Indices de Sobol pour la fonction Ishigami :
# fonctions théoriques
g0 = coeff[1]/2
g1 = sin(sort(x1))*(1 + coeff[2]/5 * pi**4)
g2 = coeff[1] * (sin(sort(x2))**2 - 1/2)
g3 = rep(0, length(x3))
g12 = rep(0, length(x1))
g13 = coeff[2] * sin(sort(x1)) * (sort(x3)**4 - 1/5 * pi**4)
g23 = rep(0, length(x3))
ish = ishigami(X, coeff)
# On centre ish
ish_centre = ish - g0
# polynomes locaux
g1_loc = predict(loess(ish_centre ~ x1, span=0.5, family="gaussian"))
g2_loc = predict(loess(ish_centre ~ x2, span=0.5, family="gaussian"))
g3_loc = predict(loess(ish_centre ~ x3, span=0.5, family="gaussian"))
options(repr.plot.width = 14, repr.plot.height = 5)
par(mfrow=c(1,3))
plot(x1, ish_centre, type = 'p')
# Prediction
lines(sort(x1), g1_loc[order(x1)], lwd=6, col="blue")
# Theorique
lines(sort(x1), g1, lwd=3, col="red")
legend(x='topleft', legend=c("Estimé", "Theorique"), col=c("blue","red"), lwd=c(2,1),lty=1,text.font=1.5,bty='n',cex=1.4)
plot(x2, ish_centre, type = 'p')
lines(sort(x2), g2_loc[order(x2)], lwd=6, col="blue")
lines(sort(x2), g2, lwd=3, col="red")
legend(x='top', legend=c("Estimé", "Theorique"), col=c("blue","red"), lwd=c(2,1),lty=1,text.font=1.5,bty='n',cex=1.4)
plot(x3, ish_centre, type = 'p')
lines(sort(x3), g3_loc[order(x3)], lwd=6, col="blue")
lines(sort(x3), g3, lwd=3, col="red")
legend(x='topleft', legend=c("Estimé", "Theorique"), col=c("blue","red"), lwd=c(2,1),lty=1,text.font=1.5,bty='n',cex=1.4)
On se rend compte qu'on ne peut pas si facilement reproduire la méthode précédente dans les cas d'intéractions du second ordre. On va donc utiliser une méthode un peu plus "manuelle". Par exemple pour simuler g13 :
Vu que g13 est la seule fonction non nulle, on ne représentera que celle-ci.
g13_chap = c()
x1 = c()
x3 = c()
nb_point = 300
#Initialisation des valeurs fixées prédéfinies
X1_grid <- X3_grid <- seq(from = -pi,to=pi,length.out =nb_point)
#Voir avancement du programme
c=0
pb <- txtProgressBar(min = 1, max = length(X1_grid), style = 3)
for (i in X1_grid) {
c=c+1
Sys.sleep(0.001)
setTxtProgressBar(pb, c)
#on fixe X1
X1 = rep(i,nb_point)
for (j in X3_grid) {
#On fixe X3
X3 = rep(j,nb_point)
#on définit X2 qu'on fait varier
X2 = runif(nb_point, min=-pi, max=pi)
contenu = c(X1,X2,X3)
#on remplit la matrice pour X1 et X3 fixés
X_fixe = matrix(contenu, nrow = nb_point, ncol=3)
Ish_X_fixe = ishigami(X_fixe, coeff)
g1_aux = sin(sort(X1))*(1 + coeff[2]/5 * pi**4)
g2_aux = coeff[1] * (sin(sort(X2))**2 - 1/2)
g13_chap_ij = mean(Ish_X_fixe)-g0-g1_aux[1]-g2_aux[1]
g13_chap = c(g13_chap, g13_chap_ij)
#pour l'affichage
x1 = c(x1, i)
x3 = c(x3, j)
}
}
|======================================================================| 100%
defaultW <- getOption("warn")
options(warn = -1)
fig <- plot_ly(x = ~X[,1], y = ~X[,3], z = ~ish_centre, mode='markers', marker = list(symbol = 'circle', sizemode = 'diameter', size = 2), type='scatter3d')
fig <- fig %>% add_trace(x = ~x1,
y = ~x3,
z = ~g13_chap, type = 'mesh3d')
fig
options(warn = defaultW)
Si l'image ci-dessus ne s'affiche pas, cf l'image .png envoyé dans le mail
Les propriétés attendues d'un plan d'expériences numériques sont :
Nous allons calculer la discrepance pour différents plan d'expériences.
critere=0
for (i in 1:10000)
{
x = matrix(runif(n=2*9),nrow=9,ncol=2) # Plan pur MonteCarlo
if (mindist(x) > critere) # selection du plan courant qui maximise
# la distance minimale entre 2 de ses points
{
x1 = x # Plan Maximin
critere = mindist(x)
}
}
dis_MC = discrepancyCriteria(x,type=c('M2','C2'))
cat("Critères discrépance M2 et C2 pour le plan pur Monte Carlo = ",
dis_MC$DisM2, " - ", dis_MC$DisC2,"\n")
dis_Maximin = discrepancyCriteria(x1,type=c('M2','C2'))
cat("Critères discrépance M2 et C2 pour le plan Maximin = ",
dis_Maximin$DisM2, " - ", dis_Maximin$DisC2,"\n")
x2 = factDesign(3,3)
dis_Fact = discrepancyCriteria(x2$design,type=c('M2','C2'))
cat("Critères discrépance M2 et C2 pour le plan Factoriel = ",
dis_Fact $DisM2, " - ", dis_Fact $DisC2,"\n")
# Création d'une suite de Sobol de 9 points en dimension 3
x3 = sobol(9,3)
dis_Sobol = discrepancyCriteria(x3,type=c('M2','C2'))
cat("Critères discrépance M2 et C2 pour la suite de Sobol = ",
dis_Sobol$DisM2, " - ", dis_Sobol$DisC2,"\n")
Critères discrépance M2 et C2 pour le plan pur Monte Carlo = 0.2656244 - 0.2316414 Critères discrépance M2 et C2 pour le plan Maximin = 0.1159467 - 0.1103984 Critères discrépance M2 et C2 pour le plan Factoriel = 0.3828553 - 0.3167492 Critères discrépance M2 et C2 pour la suite de Sobol = 0.1644181 - 0.1305912
Le plan à discrepance faible avec une suite de Sobol nous permet de minimiser le critère de discrepance. De plus, ce plan permet la séquentialité car il n'est pas nécessaire de recalculer le plan total lorsque l'on rajoute des nouveaux points. C'est donc ce plan que nous utiliserions.
# Nombre d'évaluations de f
N = 100
pex = (optimumLHS(N, 3)-0.5)*2*pi
options(repr.plot.width = 5, repr.plot.height = 5)
plot(pex, main="LHS")
# Evaluation f selon le plan pex
pex <- data.frame(pex)
y <- ishigami(pex, coeff)
# Modèle de Krigeage
m.krigeage <- km(design=pex, response=y,control=list(trace=FALSE))
options(repr.plot.width = 5, repr.plot.height = 10)
plot(m.krigeage)
#Définition de la fonction qui calcule le Q2 :
q2 <- function(y, y_pred){
return ( 1- ( (sum((y-y_pred)^2)) / sum((mean(y)-y_pred)^2) ) )
}
cross_validation <- function(pex,nb_fold){
# Evaluation de la capacité de prédiction par validation croisée
#Randomly shuffle the data
pex_sample<-pex[sample(nrow(pex)),]
#Create 10 equally size folds
folds <- cut(seq(1,nrow(pex_sample)),breaks=nb_fold,labels=FALSE)
#Perform k fold cross validation
xval_err = 0
Q2 = 0
vec_y_test = c()
vec_y_test_pred = c()
for(i in 1:nb_fold){
#Séparation des données pour la cross validation
testIndexes <- which(folds==i,arr.ind=TRUE)
testData <- pex_sample[testIndexes, ]
trainData <- pex_sample[-testIndexes, ]
#Calcul du modèle sur le jeu train
y_train = ishigami(trainData, coeff)
m.krigeage <- km(design=trainData, response=y_train,control=list(trace=FALSE))
#prediction sur le jeu de test
y_test_pred = predict.km(m.krigeage, testData, "UK", se.compute=FALSE,
checkNames=FALSE)$mean
vec_y_test_pred = c(vec_y_test_pred,y_test_pred)
y_test = ishigami(testData, coeff)
vec_y_test = c(vec_y_test,y_test)
e = y_test_pred-y_test
xval_err = xval_err + e%*%e/length(y_test)
Q2 = Q2 + q2(y_test, y_test_pred)
}
rmse = sqrt(xval_err/nb_fold)
Q2 = Q2/nb_fold
return(list(RMSE = rmse, Q2=Q2, Vec_pred = vec_y_test_pred, Vec_reel = vec_y_test))
}
list_result = cross_validation(pex,10)
cat("Valeur de la RMSE : " , list_result$RMSE) #doit être le plus petit possible
Valeur de la RMSE : 1.304057
cat("Valeur du Q2 : " ,list_result$Q2) # doit etre le plus proche de 1
Valeur du Q2 : 0.8491658
#Graphe des valeurs prédites VS valeurs réelles :
options(repr.plot.width = 5, repr.plot.height = 5)
y_test = list_result$Vec_reel
y_test_pred = list_result$Vec_pred
plot(y_test, y_test_pred, main = "Valeurs prédites vs Valeurs réelles", xlab="Reelles", ylab="Predites")
reg <- lm(y_test_pred ~ y_test)
abline(reg, col='red')
lines(y_test,y_test, col='blue')
legend(x='bottom', legend=c("Reg Predite", "Reg Reelle"), col = c('red','blue'), cex= 0.9, lty=1, bty='n')
m.krigeage <- km(design=pex, response=y,control=list(trace=FALSE))
kriging.mean <- function(Xnew, m){
predict.km(m, Xnew, "UK", se.compute = FALSE, checkNames = FALSE, control=list(trace=FALSE))$mean
}
SA.metamodel <- fast99(model = kriging.mean, factors = 3, n = 1000,
q = "qunif", q.arg = list(min = -pi, max = pi), m = m.krigeage)
#SA_test <- fast99(model = ishigami.fun, factors = 3, n = 1000,
# q = "qunif", q.arg = list(min = -pi, max = pi))
options(repr.plot.width = 5, repr.plot.height = 5)
plot(SA.metamodel, cex=0.8)
Interprétation :
La variable qui a le plus d’influence sur la variance de la sortie (c'est à dire avec l'indice total le plus élevé) est la variable X1, avec un indice total d'environ 60% et un indice du premier ordre d'environ 30%. (X1 explique à elle seule environ 30% de la variance de Y)
La variable X2 n’intervient que seule dans la fonction Ishigami, cela se retrouve sur le graphique ci-dessus où son indice du premier ordre est quasiment équivalent à son indice total. (environ 40% de la variance de Y est expliquée par X2)
La variable X3 n’a aucune influence seule, mais a une influence relativement importante en intéraction et notamment avec X1 (comme on le retrouve d'ailleurs dans la formule d'Ishigami).
On en déduit que la variance de Y est due pour 40% à X2, 30% à X1 et 30% à l’interaction entre X1 et X3.
cat("Indices de Sobol pour la fonction Ishigami : \n")
(ishigami.sobol.indices(coeff))
print(SA.metamodel)
Indices de Sobol pour la fonction Ishigami :
Call: fast99(model = kriging.mean, factors = 3, n = 1000, q = "qunif", q.arg = list(min = -pi, max = pi), m = m.krigeage) Model runs: 3000 Estimations of the indices: first order total order X1 0.296888140 0.5186778 X2 0.455851334 0.5031961 X3 0.004771931 0.2350741
coeff_ordre1_metam=SA.metamodel$D1/SA.metamodel$V
coeff_tot_metam = 1-((SA.metamodel$Dt)/SA.metamodel$V)
coeff_ordre1=ishigami.sobol.indices(coeff)$Coeff_Ordre1
coeff_ordre_tot=ishigami.sobol.indices(coeff)$Coeff_Ordre_Tot
cat("Différence absolue les valeurs estimées avec le métamodèle et les valeurs théoriques des coefficients d'ordre 1 :\n")
cat(abs(coeff_ordre1_metam-coeff_ordre1), '\n \n')
cat("Différence les valeurs estimées avec le métamodèle et les valeurs théoriques des coefficients d'ordre total :\n")
cat(abs(coeff_tot_metam-coeff_ordre_tot))
Différence absolue les valeurs estimées avec le métamodèle et les valeurs théoriques des coefficients d'ordre 1 : 0.01701705 0.01344019 0.004771931 Différence les valeurs estimées avec le métamodèle et les valeurs théoriques des coefficients d'ordre total : 0.03891105 0.06078491 0.008609577
n = 4000
X = matrix(NA, n, 3) #Creation d'une matrice vide
x1 = runif(n,-pi,pi)
x2 = runif(n,-pi,pi)
x3 = runif(n,-pi,pi)
X[,1] = x1
X[,2] = x2
X[,3] = x3
#initialisation des coefficients A et B de la fonciton d'Ishigami:
coeff = c(7, 0.1)
ish_krig = kriging.mean(X,m.krigeage)
# On centre ish_krig
ish_krig_centre = ish_krig - g0
# polynomes locaux
g1_loc_krig = predict(loess(ish_krig_centre ~ x1, span=0.5, family="gaussian"))
g2_loc_krig = predict(loess(ish_krig_centre ~ x2, span=0.5, family="gaussian"))
g3_loc_krig = predict(loess(ish_krig_centre ~ x3, span=0.5, family="gaussian"))
options(repr.plot.width = 14, repr.plot.height = 5)
par(mfrow=c(1,3))
plot(x1, ish_krig_centre, type = 'p')
# Prediction
lines(sort(x1), g1_loc_krig[order(x1)], lwd=6, col="blue")
# Theorique
lines(sort(x1), g1, lwd=3, col="red")
legend(x='topleft', legend=c("Estimé (krigeage)", "Theorique"), col=c("blue","red"), lwd=c(2,1),lty=1,text.font=1.5,bty='n',cex=1.4)
plot(x2, ish_krig_centre, type = 'p')
lines(sort(x2), g2_loc_krig[order(x2)], lwd=6, col="blue")
lines(sort(x2), g2, lwd=3, col="red")
legend(x='top', legend=c("Estimé (krigeage)", "Theorique"), col=c("blue","red"), lwd=c(2,1),lty=1,text.font=1.5,bty='n',cex=1.4)
plot(x3, ish_krig_centre, type = 'p')
lines(sort(x3), g3_loc_krig[order(x3)], lwd=6, col="blue")
lines(sort(x3), g3, lwd=3, col="red")
legend(x='topleft', legend=c("Estimé (krigeage)", "Theorique"), col=c("blue","red"), lwd=c(2,1),lty=1,text.font=1.5,bty='n',cex=1.4)
Avec seulement 100 simulations, on a des bonnes approximations. Cela montre l'intérêt de passer par des métamodèles.
# Nombre d'évaluations de f
N_seq = c(20,50,100,150)
nb_eval = 5
nb_fold = 10
df_rmse <- list()
df_Q2 <- list()
i= 1
for (Ni in N_seq) {
vec_rmse = c()
vec_Q2 = c()
for (k in 1:nb_eval){
pex_i = (optimumLHS(Ni, 3)-0.5)*2*pi
# Evaluation f selon le plan pex
pex_i <- data.frame(pex_i)
list_result_i = cross_validation(pex_i,10)
vec_rmse <- c(vec_rmse, list_result_i$RMSE[1])
vec_Q2 <- c(vec_Q2, list_result_i$Q2)
}
print(vec_rmse)
df_rmse[[i]] <- vec_rmse
df_Q2[[i]] <- vec_Q2
i = i+1
}
#options(repr.plot.width = 5, repr.plot.height = 5)
#plot(pex, main="LHS")
[1] 4.594274 4.969337 3.978402 3.344505 3.276791 [1] 2.950537 3.295803 2.225401 2.274905 3.703246 [1] 1.375695 1.154038 1.461013 1.339915 1.335368 [1] 0.5825002 0.8992492 0.8225002 0.6376081 0.7362206
boxplot(df_rmse,main='Boxplot de RMSE selon la taille du LHS',
names=c(20,50,100,150),
xlab='Taille LHS',ylab='RMSE',col='lightblue')
Plus le nombre de simulations est élevé, plus la RMSE diminue.